14. Análisis de Asociación de Variables Aleatorias#
En esta sección, se presentan tests para analizar la asociación entre variables aleatorias continuas (ANOVA), o la asociación entre variables aleatorias discretas o categóricas (Chi-cuadrado, Bowker).
14.1. ¿Por qué ANOVA?#
En los test de hipótesis estudiados en sesiones anteriores, consideramos el uso de la distribución \(t\)- student para analizar las medias muestrales. En lo que sigue utilizaremos el Análisis de la Varianza (ANOVA) o también llamado análisis de factores, para estudiar el efecto de uno o más factores (cada uno con dos o más niveles) sobre la media de una variable continua.
Test de medias para dos poblaciones normales con la misma varianza desconocida (caso común)
Sean \(X_1,\cdots,X_n\) e \(Y_1,\cdots,Y_m\) muestras independientes de poblaciones normales con medias desconocidas \(\mu_x\) y \(\mu_y\) y misma varianza desconocida \(\sigma^2\). Consideremos el test de hipótesis:
del Corolario del Teo de Fisher-Cochran se cumple:
donde
de manera que se rechaza \(H_0\) si
En el caso en que se requiera comparar más de 2 grupos, o examinar el efecto de 1, 2 o mas factores, este procedimiento se vuelve ineficiente, porque no queremos hacer un montón de t-tests para cada par.
Además, family-wise error rate (la tasa de error de la familia, error global) que es la probabilidad de cometer al menos un error de Tipo I en múltiples pruebas estadísticas realizadas en los mismos datos aumenta. Para \(c\) tests se calcula como:
con c=2 tests, el error de tipo I es 0.0975 (alrededor de 2*0.05=0.1)
con c=3 tests, el error de tipo I es 0.143 (alrededor de 3*0.05=0.15)
con c=10 tests, el error de tipo I es 0.40 (alrededor de 10*0.05=0.5)
Por lo tanto, es mejor abordar con el modelo del Análisis de la Varianza (ANOVA).
14.2. Repaso#
14.2.1. La distribución chi-cuadrado#
Sean \(Z_1,\cdots, Z_n\, v.a. i.i.d. \, \sim {\it N}(0,1)\) entonces
donde \(n\) son los grados de libertad de la distribución.
Propiedades de la distribución \(\chi^2\):
(i) Propiedad aditiva: si \(X_1\) y \(X_2\) son dos v.a. independientes distribuidas \(\chi^2\) de \(n_1\) y \(n_2\) grados de libertad, entonces
(ii) Esperanza y Varianza:
suppressMessages(library(dplyr))
suppressMessages(library(plotly))
suppressMessages(library(ggplot2))
suppressMessages(library(rmarkdown))
vec <- seq(0,100,by=0.1)
params <- c(1:5, seq(10, 100, by=10))
pvec <- list()
for (i in 1:length(params)){
pvec[[i]] <- dchisq(vec,df=params[i],ncp=0)
}
steps <- list()
fig <- plot_ly(width=600,height=600) %>% layout(title = "\n \n Densidad de Probabilidad Chi-cuadrado", yaxis = list(range=c(0,0.3)))
for (i in 1:length(params)){
fig <- add_lines(fig, x=vec, y=pvec[[i]], visible=if (i==1) TRUE else FALSE, mode='lines', line=list(color='blue'), showlegend=FALSE)
steps[[i]] = list(args = list('visible', rep(FALSE, length(params))), label=params[i], method='restyle')
steps[[i]]$args[[2]][i] = TRUE
}
fig <- fig %>% layout(sliders = list(list(active=0, currentvalue = list(prefix = "df: "), steps=steps)))
fig
14.2.2. La distribución F#
Sean \(X \sim \chi_n^2\) e \(Y \sim \chi_m^2\) dos v.a. independientes \(\chi^2\) de grados de libertad \(n\) y \(m\) respectivamente, entoncese se define:
donde \(F_{n,m}\) es la distribución \(F\) de \(n\) y \(m\) grados de libertad. También se nota \(F(n,m)\).
vec <- seq(0,5,by=0.01)
params <- seq(1,20,by=1)
pvec <- list()
for (i in 1:length(params))
for (j in 1:length(params)){
k = length(params)*(i-1) + j
pvec[[k]] <- df(vec,df1=params[i],df2=params[j],ncp=0)
}
steps1 <- list()
steps2 <- list()
fig <- plot_ly(width=600,height=600) %>% layout(title = "\n \n Densidad de Probabilidad F(n, m)",
yaxis = list(range=c(0,1)))
for (i in 1:length(params)){
for (j in 1:length(params)){
k = length(params)*(i-1) + j
fig <- add_lines(fig, x=vec, y=pvec[[k]],
visible=if ((i==1) && (j==1)) TRUE else FALSE,
mode='lines', line=list(color='blue'), showlegend=FALSE)
steps2[[j]] = list(args = list('visible', rep(FALSE, length(params)*length(params))),
label=params[j], method='restyle')
steps2[[j]]$args[[2]][k] = TRUE
steps1[[i]] = list(args = list('visible', rep(FALSE, length(params)*length(params))),
label=params[i], method='restyle')
steps1[[i]]$args[[2]][k] = TRUE
}
}
fig <- fig %>% layout(sliders =
list( list(active=0, currentvalue = list(prefix = "df1 (n): "), pad = list(t=20), steps=steps1),
list(active=0, currentvalue = list(prefix = "df2 (m): "), pad = list(t=100), steps=steps2)))
fig
Percentil
Sea \(F_{\alpha,n,m}\) tal que \(P\{F > F_{\alpha,n,m} \} = \alpha\), \(F_{\alpha,n,m}\) es el percentil \(100(1-\alpha)\) de la distribución \(F_{n,m}\).
Relación entre las distribuciones F y t
Si \(X \sim t_{n}\) entonces \(X^2 \sim F_{1,n}\). Por qué?
Por construcción, si \(X \sim t_{n}\), entonces
Entonces
14.2.3. La distribución de la media y varianza muestral del caso Normal#
Teorema de Fisher-Cochran
Sean \(X_1,\cdots,X_n\) v.a. i.i.d. \({\cal N}(\mu,\sigma^2)\) entonces la media y varianza muestral cumplen:
\(\begin{equation} \begin{array}{lcll} (i) & \bar{X} &\sim& {\cal N}(\mu, \frac{\sigma^2}{n})\\ \\ (ii) & {\displaystyle \frac{(n-1)S^2}{\sigma^2}}& \sim& \chi_{n-1}^2 \\ \\ (iii)& \bar{X} &{\mathrel \perp} & S^2 \quad \text{(independentes)}\\ \end{array} \end{equation}\)
\(S^2\) es la varianza muestral:
Por qué en \( {\displaystyle \frac{(n-1)S^2}{\sigma^2}} \sim \chi_{n-1}^2\), el grado de liberdad es \(n-1\)?
14.3. ANOVA unidireccional (ANOVA de un factor)#
Suponga que tenemos \(m\) tratamientos distintos, donde el resultado de aplicar el tratamiento \(i\) a un individuo es una v.a. normal. Estamos interesados en probar la hipótesis de que todos los tratamientos tienen el mismo efecto. Para ello se aplica cada tratamiento a una muestra distinta de \(n\) individuos y se analizan los resultados.
Formalmente, sean \(m\) muestras (grupos) independientes de tamaño \(n\):
donde \(\mu_i\) y \(\sigma^2\) son desconocidos. Estamos interesados en probar:
14.3.1. Supuestos (y robustez)#
Consideramos sus supuestos y robustez aquí:
Normalidad: cada muestra (o v.a. \(X_{ij}\)) debe provenir de una distribución normal. Utilizar por ejemplo, test de Shapiro-Wilk. Utilizar Q-Q plot para descartar la presencia de outliers.
El ANOVA unidireccional se considera una prueba robusta contra el supuesto de normalidad, cuando el tamaño de cada muestra (o grupo) no es pequeño (y los tamaños de grupos son iguales); si no, consideremos test nonparamétrico como Kruskal-Wallis H Test.
Varianza común (homogeneidad): misma varianza en los grupos. Se puede probar la igualdad de varianzas con el test de Levene o con el test de Brown-Forsythe.
Si este supuesto falla, consideraremos test nonparamétrico como Kruskal-Wallis H Test.
Observaciones independientes. Las poblaciones de los grupos son independientes y las observaciones de cada grupo son independientes. Este supuesto no puede ser probado con ningún estadístico, es una consideración de diseño
La violación de esto es muy grave (por ejemplo, la tasa de error de Tipo I está sustancialmente inflada).
Idea general de la derivación ANOVA
Se trata de construir dos estimadores de la varianza común, el primer estimador no depende de si \(H_0\) es cierto o no. En cambio, el segundo estimador asume que \(H_0\) es cierto, en caso contrario este estimador sobreestima a \(\sigma^2\). El test compara ambos estimadores y rechaza \(H_0\) cuando la tasa entre el segundo y primer estimador es suficientemente grande.
14.3.2. El primer estimador de \(\sigma^2\) que no depende de \(H_0\) (\(SS_W\))#
Definimos:
Recordemos la propiedad:
Así que:
Denominaremos
(\(SS_W\) es “within samples sum of squares”, o “residual sum of squares”)
Entonces
y luego
con lo cual tenemos nuestro primer estimador insesgado (por qué insesgado?) de \(\sigma^2\):
que se conoce como “error cuadrático medio” (mean square error), “cuadrado de la media de los residuos” (residual mean square) o “varianza intragrupos” (within-group variance). Notar que no considera ningun supuesto sobre \(H_0\).
14.3.3. El segundo estimador de \(\sigma^2\) que depende de \(H_0\) (\(SS_B\))#
En este caso suponemos que se cumple \(H_0\), es decir
De manera de que las medias muestrales cumplen:
entonces:
Denominaremos
(\(SS_B\) es “between samples sum of squares”, o “model sum of squares”.)
Entonces, si \(H_0\) es cierto, se cumple
y luego
con lo cual tenemos nuestro segundo estimador insesgado (por qué insesgado?) de \(\sigma^2\):
que se conoce como “cuadrado de la media entre muestras” (between samples mean square), “cuadrado de la media del modelo” (model mean square) o “varianza entregrupos” (between-group variance).
14.3.4. Test estadístico#
Se define como test estadístico (Test Statistics):
Cuándo rechazamos la hipótesis nula?
Se rechazará la hipótesis nula cuando TS sea suficientemente grande.
¿Cuál es la distribución de TS?
Demostración:
Sea entonces \(F\sim F_{m-1,nm-m}\) y \(F_{m-1,nm-m,\alpha}\) tal que:
Se calcula entonces el valor de \(TS\), si \(TS > F_{m-1,nm-m,\alpha}\) se rechaza \(H_0\).
14.3.5. Identidad de suma de cuadrados#
Definimos \(SS_{total} = SS_{T}\) (total sum of squares) como:
Cumple:
Así que tenemos:
14.3.6. Ejemplo (industria de automóviles)#
Una industria de automóviles desea probar la calidad de 3 tipos de combustibles. Para ello observa el rendimiento (millas por 10 galones de combustible) de los 3 tipos de combustibles en 15 vehículos idénticos (5 con cada tipo de combustible). (Asumimos que los supuestos de ANOVA se mantienen aquí, pero hay que examinar los supuestos en tu análisis de los datos.)
Cuál es la hipótesis nula aquí?
Veamos primero la distribución de los datos:
library(tidyr)
library(ggplot2)
library(dplyr)
library(cowplot)
options(repr.plot.width=10, repr.plot.height=5)
options(digits=4)
gas1 <- c(220,251,226,246,260)
gas2 <- c(244,235,232,242,225)
gas3 <- c(252,272,250,238,256)
df <- do.call(rbind, Map(data.frame, gas1=gas1, gas2=gas2, gas3=gas3))
df <- df %>% pivot_longer(cols=c('gas1', 'gas2', 'gas3'),names_to='combustible', values_to='milla')
df %>% group_by(combustible) %>% summarise(M=mean(milla), SD=sd(milla))
#mean_sdl: computes the mean plus or minus a constant times the standard deviation. By default mult=2
p1 <- ggplot(df, aes(x=combustible, y=milla)) +
geom_violin(trim=FALSE) +
geom_dotplot(binaxis='y', stackdir='center', dotsize=1, binwidth=1.5) +
stat_summary(fun.data=mean_sdl, fun.args=list(mult=1), geom="pointrange", color="red", lwd=1) +
theme(text=element_text(size=14))
p2 <- ggplot(df, aes(x=combustible, y=milla)) + geom_boxplot() + theme(text=element_text(size=14))
plot_grid(p1, p2, align="h", nrow=1, ncol=2, rel_widths= c(1/2, 1/2))
| combustible | M | SD |
|---|---|---|
| <chr> | <dbl> | <dbl> |
| gas1 | 240.6 | 16.965 |
| gas2 | 235.6 | 7.701 |
| gas3 | 253.6 | 12.280 |
En tu tarea, tienes que desarrollar el test. Aquí, mostramos el resultado de R:
fm <- aov(milla ~ combustible, data=df)
summary(fm)
Df Sum Sq Mean Sq F value Pr(>F)
combustible 2 863 432 2.6 0.12
Residuals 12 1992 166
14.3.7. Caso de muestras con diferentes tamaños#
¿Qué ocurre si las muestras no son del mismo tamaño?
donde \(N = \sum_{i=1}^m n_i - m\). Y entonces,
por otra parte
y se definen
Finalmente, si \(TS > F_{\alpha,m-1,N}\) se rechaza \(H_0\).
14.4. Métodos de comparaciones múltiples#
Si el test ANOVA es estadísticamente significativo, lo único que podemos concluir es que una o mas de las diferencias entre grupos es significativa, pero no sabemos que grupos son los que difirien. Para ello es necesario realizar “Comparaciones múltiples”. Se trata de comparar cada grupo con otro grupo, o un promedio de grupos se puede comparar con otros. Se consideran dos posibilidades:
Comparaciones planeadas: existe interés por algunos grupos en particular
Comparaciones post hoc: no existen hipótesis específicas de algunos grupos. Nos centramos en este caso aquí.
Aquí presentamos algunos métodos de comparaciones múltiples que son los más comúnmente usados.
14.4.1. Bonferroni#
Utiliza el mismo estadístico que LSD, pero calcula la significancia \(\alpha\) considerando que se realizan c comparaciones.
Dos posibilidades:
Se calcula un nuevo valor de \(\alpha\) para mantener el error de tipo I global:
Por ejemplo, si \(m=4, c=6\) y \(\alpha=0.05\), entonces \(\alpha_{nuevo} = 0.0083\)
Se corrige el p-value multiplicando por el número de comparaciones, y lo compara con el \(\alpha\) original (e.g., 0.05):
Bonferroni puede ser demasiado conservador (\(p_{Bonferroni}\) bastante grande, o \(\alpha_{nuevo}\) bastante pequeño), lo que dificulta detectar diferencias que realmente existen (reduciendo su potencia). Pero puede garantizar el error de Tipo I, y es facíl entender la corrección. También lo usamos en comparaciones planeadas.
14.4.2. Tukey y Tukey-Kramer#
Tukey: mismo tamaño muestral n
Se define el siguiente estadístico:
Tukey-Kramer: distintos tamaños muestrales
donde
En ambos casos \(q_N\) sigue una distribución de rangos student (studendized range distribution) que es similar a t-student (studentized t-test) pero diferente. Usabmos el \(\alpha\) original (e.g., 0.05). Segimos el proceso de de valor crítico o p-value.
El Tukey es menos conservador que Bonferroni. Tiene un control estricto sobre error de Tipo I (menos estricto que Bonferroni). Tiene buena potencia.
Algunas preguntas
Si el test ANOVA es estadísticamente significativo, implica que ¿existe al menos un grupo cuyo comportamiento difiere de los otros?
No necesariamente, podría ocurrir que sea la combinación de grupos que provoque las diferencias. El test de Scheffe permite pesquisar las diferencias subyacentes.
¿Es posible encontrar diferencias siginificativas con las comparaciones múltiples, aunque ANOVA no haya sido significativo?
Si, es posible, debido a que los tests de comparaciones múltiples están mas enfocados (tienen más potencia) que ANOVA.
¿Es útil el resultado del test ANOVA?
Claro, cuando corresponde a la hipótesis en estudio!
A menudo, sus preguntas experimentales están más enfocadas en algunos grupos. En estos casos, puede saltar directamente a los tests de comparaciones múltiples (aunque algunas personas no están de acuerdo con esta afirmación).
Pero tenga en cuenta que todos los cálculos de comparación múltiple utilizan MSE de la tabla ANOVA.
14.5. Análisis de asociación en v.a. discretas (o categóricas)#
En esta sección se presentan algunos tests que permiten analizar la asociación entre variables aleatorias discretas.
14.5.1. Tabla de contingencia#
Consideremos que tenemos observaciones de dos variables discretas o categóricas (o bien variables continuas agrupadas en clases) \(X\) y \(Y\) referidas a una muestra. La tabla de contingencia contiene los efectivos correspondientes a cada pareja de valores (o clases) de ambas variables:
Cada linea y cada columna corresponden a una submuestra particular. La fila de índice \(h\) es la distribución sobre \( d_1,\ldots,d_s\) de los individuos para los cuales la variable \(X\) toma el valor \(c_h\). La columna de índice \(k\) es la distribución sobre \(c_1,\ldots,c_r\) de los individuos para los cuales las variable \( Y\) toma el valor \(d_k\).
Los totales de renglones y columnas en la tabla \(n_{h\bullet}, n_{\bullet k}\) (etc.) se denominan frecuencias marginales.
Los valores en las celdas internas \(n_{hk}\) (etc.) se denominan frecuencias observadas.
14.5.2. Test de Chi-cuadrado de independencia para tablas de contingencia#
Tenemos observaciones de dos variables discretas o categóricas (o bien variables continuas agrupadas en clases) \(X\) y \(Y\) referidas a una muestra de distribuciones desconocidas. Las observaciones de \(X\) son independendes, y las observaciones de \(Y\) son independendes. Nos interesa si \(X\) y \(Y\) son independientes. Tenemos las hipótesis:
donde
\(P_{hk}\) se refiere a la probabilidad conjunta de que \(X\) sea \(c_h\) y \(Y\) sea \(d_k\)
\(p_{h}\) se refiere a la probabilidad marginal de que \(X\) sea \(c_h\)
\(q_{k}\) se refiere a la probabilidad marginal de que \(Y\) sea \(d_k\)
Ya que no conocemos sus distribuciones, estimamos \(P_{hk}, p_{h}, q_{k}\) a partir de la muestra, utilizando la tabla de contingencia:
(Por qué podemos estimar así?)
Asi que bajo \(H_0\), queremos que los valores de \(\hat P_{hk}\) y \(\hat p_{h} \hat q_{k}\) estén cercas:
Recuerda que \(n_{hk}\) es la frecuencia observada (\(O_{hk}\)), y aquí el producto \(\displaystyle \frac{n_{h \bullet} n_{\bullet k}}{n}\) se denomina la frecuencia esperada (\(E_{hk}\)). Definimos un estadístico (v.a.) que representa la distancia entre ellos para todas combinaciones de valores:
Tenemos un teorema que dice bajo \(H_0\), este estadístico \(TS\) converge hacia una distribución \(\chi^2_{(r-1)(s-1)}\) cuando \( n\) tiende al infinito.
Advertencia
Para que la aproximación chi-cuadrado sea válida, cada frecuencia esperada debe ser de al menos 5 (esta prueba no es válida para muestras pequeñas). Si alguna frecuencia esperada sea menor que 5, consideramos combinar algunas clases o usamos la prueba exacta de Fisher / Fisher’s Exact Test o Monte Carlo simulación (chisq.test(table, simulate.p.value=T)).
Por qué el grado de libertad es \((r-1)(s-1)\)?
Ejemplo (campus y tiempo libre)
Consideremos dos variables categóricas referidas a una muestra de los estudiantes de la UACh, una de las cuales indica el campus, y la otra indica que hace en su tiempo libre. Se trata de analizar si existe una dependencia entre el campus y actividades en su tiempo libre. Supongamos que la tabla de contingencia observada es la siguiente (los números son inventados):
En tu tarea, tienes que desarrollar el test. Aquí, mostramos el resultado de R:
#usando el test predifinido en R
table <- rbind(c(68, 56, 32), c(52, 72, 20))
chisq.test(table, correct=FALSE)
Pearson's Chi-squared test
data: table
X-squared = 6.4, df = 2, p-value = 0.04
14.5.3. Test de Bowker (o McNemar-Bowker)#
Este test se utiliza para analizar la simetría en una tabla de contingencia construida a partir de dos muestras dependientes \(X, Y\) (donde las observaciones de cada muestra son independentes), típicamente pre y post un tratamiento.
Las hípotesis de este test es:
donde
\(p_{ij}\) se refiere a la probabilidad de que \(X\) sea \(i\) y \(Y\) sea \(j\) (usamos \(i\) o \(j\) para referir a un valor de una clase o categoría, aunque strictamente debemos escribir \(c_i\) o \(c_j\))
\(p_{ji}\) se refiere a la probabilidad de que \(X\) sea \(j\) y \(Y\) sea \(i\)
El estadístico del test de Bowker se construye como:
De manera que bajo la hipótesis nula de simetría en la tabla de contingencia, es decir que no hay efecto del tratamiento, se cumple asintóticamente que:
Si se rechaza la hipótesis nula (cuándo?), y además se observan más efectivos sobre la diagonal (upper off-diagonal) que bajo la diagonal (lower off-diagonal) se puede concluir que hay un efecto positivo del tratamiento en el criterio en análisis.
Advertencia
Para que la aproximación chi-cuadrado sea válida, se recomienda que \(n_{ij} + n_{ji} \geq 10, \forall i, j, i \neq j\) (o \( n_{ij} + n_{ji} \geq 25\) según otro fuente). Si no se cumple, consideramos combinar algunas clases o usamos el mid-P McNemar (binomial) test.
Si \(n_{ij} + n_{ji} = 0\) en algún par, lo ignoramos en la calculación de B.
Ejemplo
Considere los datos de 170 estudiantes sometidos a una evaluación sobre sus habilidad de cohesión en la redacción de textos, antes y después de una intervención didáctica. La evaluación utiliza una escala de 4 niveles desde no logrado a totalmente logrado. Se pide evaluar si la intervención ha tenido efecto o no en la habilidad estudiada.
#fila (i): pre-evaluación
#columna (j): post-evaluación
tabla <- rbind(c(0,6,9,6), c(6,6,9,21), c(9,12,18,55), c(4,0,3,6))
print(tabla)
B = 0
efsd = efbd = 0
for (i in 1:3){
for (j in (i+1):4){
B = B + (tabla[i,j] - tabla[j,i])^2 / (tabla[i,j] + tabla[j,i])
efsd = efsd + tabla[i,j]
efbd = efbd + tabla[j,i]
}
}
cat("B:", B, "\n")
cat("sobre diag vs. bajo diag:", c(efsd, efbd), "\n")
k = 4*3/2
q = qchisq(0.95, k)
print(q)
p <- 1-pchisq(B, k)
print(p)
[,1] [,2] [,3] [,4]
[1,] 0 6 9 6
[2,] 6 6 9 21
[3,] 9 12 18 55
[4,] 4 0 3 6
B: 68.44926
sobre diag vs. bajo diag: 106 34
[1] 12.59159
[1] 8.500978e-13
Se rechaza la hipótesis nula de simetría en favor de un efecto positivo de la intervención didáctica.
##usando el test predefinido en R
mcnemar.test(tabla)
McNemar's Chi-squared test
data: tabla
McNemar's chi-squared = 68.449, df = 6, p-value = 8.501e-13